Diffusion limited aggregation as a Markovian process. 
Part I: bond-sticking conditions 
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Cylindrical lattice Diffusion Limited Aggregation (DLA), with a narrow width N, is solved using 
a Markovian matrix method. This matrix contains the probabilities that the front moves from one 
configuration to another at each growth step, calculated exactly by solving the Laplace equation 
and using the proper normalization. The method is applied for a series of approximations, which 
include only a finite number of rows near the front. The matrix is then used to find the weights of 
the steady state growing configurations and the rate of approaching this steady state stage. The 
former are then used to find the average upward growth probability, the average steady-state density 
and the fractal dimensionality of the aggregate, which is extrapolated to a value near 1.64. 



I. INTRODUCTION 

Diffusion limited aggregation (DLA)i has been the 
subject of extensive study since it was first introduced. 
This model exhibits a growth process that produces 
highly ramifipjd self similar patterns, which are believed 
to be fractalsH. It seems that DLA captures the essential 
mechanism in many natural growth processes, such as 
viscous fingeringcl, dielectric breakdownel, etc. It is now 
understood that the Laplace equation, which is common 
to all of these processes and to DLA, has a major role 
in the resemblance between them. One of the interest- 
ing features of DLA is that there are no parameters to 
fine-tune in order to obtain a fractal. It thus differs from 
ordinary critical phenomena, and belongs to the class of 
self organized criticality (SOC)El. In spite of the appar- 
ent simplicity of the model, an analytic solution is still 
unavailable. Particularly, the exact value of the fractal 
dimension is not known. 

In DLA there is a seed cluster of particles fixed some- 
where. A particle is released at a distance from the clus- 
ter, and performs a random walk until it attempts to 
penetrate the fixed cluster, in which case it sticks. Then 
the next particle is released and so on. There are two 
common types of sticking conditions. The sticking con- 
dition described above is called "bond-DLA" , because it 
occurs when a particle goes into a perimeter bond. In 
"site-DLA", sticking occurs as soon as the particle ar- 
rives in a perimeter site. This paper deals with bond- 
DLA, whereas part II deals with site-DLA. The large 
scale structure of DLApis not sensitive to the type of 
sticking conditions useclEnj. 

It has been shown that bond-DLA is equivalent to the 
dielectric breakdown model (DBM) with rj = llil. DBM 
is a cellular automaton that is defined on a lattice. It 
consists of the following steps: one starts with a seed 
cluster of connected sites and with a boundary surface 
far away from it. A field $, which corresponds to the 



electrostatic potential, is found by solving the discrete 
Laplace equation on a lattice, 



= 0, 



(1.1) 



with the following boundary conditions: the aggregate 
is considered to have a constant potential that is usu- 
ally set to 0, and the potential gradient on the distant 
boundary is set to 1 in some arbitrary units (some use a 
constant potential on the distant boundary instead). In 
this paper we set the distant boundary at infinity, and ig- 
nore the exponentially small finite size corrections. After 
solving the discrete Laplace equation (1.1), the field <& 



determines the growth probabilities per perimeter bond. 
More specifically, the growth probabilities are propor- 
tional to the electric field to some power r\. The electric 
field is simply equal to the potential difference across each 
bond. Because the potential is set to on the aggregate, 
the electric-field is equal to the potential value at the 
sites lying across the perimeter bonds. Thus, 



Pb = 



l*ftl" 



£ 6 W 



(1.2) 



Here, b is the bond index. 

DLA and DBM can be grown in various geometries. 
By geometry we refer to the dimensionality of the lat- 
tice, to the shapes of the boundaries and to the details 
of the seed for growth (usually a point or a line for two 
dimensional growth). For instance, the case in which the 
distant boundary is a sphere is called radial boundary 
conditions, and the case in which the boundary is a dis- 
tant plane at the top, while the seed cluster is a parallel 
plane at the bottom, with periodic boundary conditions 
on the sides, is called cylindrical boundary conditions. 
In this paper we only consider the cylindrical case, with 
a relatively short period length (width), from N = 2 to 
about N — 7, although the method described here could 
also be used for wider cases. 
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Recently we published an exact solution to DLA in 
cylindrical geometry of width N = 2u. The present pa- 
per generalizes and extends that solution. Our approach 
follows the dynamics of the interface. The interface alone 
determines the growth probabilities at each time step, 
and whatever lies behind it is irrelevant. This is because 
the solution to the Laplace equation is unique, provided 
that the boundary conditions-are well defined. We now 
give a brief summary of Refill The characterization of 
the interface for N = 2 is simple; The interface is fully 
characterized by a single parameter (usually denoted by i 
or j), which corresponds to the height difference between 
the two columns. This height difference, referred to as 
the step size, can be infinitely large; see Fig. 0. If the 
interface is flat (J = 0), one can assume that the next 
particle will always stick on the right side, without limit- 
ing the generality of this discussion. This means that the 
step size can always be considered as nonnegative. The 
Markovian dynamics is then presented using the Master 
equation, 



Pi(t + l)=^ / E iJ P j (t), 



(1.3) 



where Pj(t) is the probability that the step size is j at 
time t, and Eij is the time independent conditional prob- 
ability that an initial step size j will become i after the 
next growth process. An example with several possible 
transitions is shown in Fig. |^. P(t) is called the state 
vector and E is called the evolution matrix. In princi- 
ple, a similar Master equation can be written for more 
complex growth situations, provided the various configu- 
rations can be indexed with a single index j. Being made 
out of conditional probabilities, the elements of the evo- 
lution matrix obey that, 



< E it j < 1, i, j = 0, ...,oo, 
>:,%/••,,. I- J=0,...,oo. 



(1.4) 



After many iterations of Eq. (1.3) the system converges 
to a fixed point P*, also called the steady state, which 
represents the asymptotic time distribution of the step 
sizes. From the steady state and the evolution matrix we 
are able to extract the average upward growth probability 
(p up )*, the average density p and the fractal dimension 
D. 

In order to obtain an analytic expression for the ele- 
ments of the evolution matrix, one must first solve the 
Laplace equation. Having found the solution s $ (m, n), 
the growth probabilities are found from Eq. (1.2). The 
denominator there, which comes from the normalization, 
is particularly simple for the special case of r\ — 1 , where 
the .discrete version of the divergence theorem implies 
thatfl 



(1.5) 



The actual growth probability into a site is then found 
from 



Psitc 



E 



Pb- 



(1.6) 



bonds into site 



The solution of the Laplace equation is now divided 
into two parts. In the first part, we solve the Laplace 
equation for the 'upper' part of space, which starts just 
above the highest particle of the aggregate and continues 
upwards to infinity. In the example of Fig. 1, this part 
contains all the rows with m > 0. As we explain below, 
this solution is completely determined by the boundary 
conditions and by the values of the potential at the row 
with m = 0, i.e. {$(0,n)}. We then solve the Laplace 
equation for the 'lower' part (to < in Fig. 1), and find 
the values of {<&(0,n)} from matching the two regimes. 
The solution in the 'upper' part is given as a combination 
of solutions of the formtl 



$(w,n) = e Km+4fen , 
with the dispersion relation 



smh(-) = ±sin(-) 



(1.7) 



(1.8) 



and with the discrete set of allowed values k\ = ■jj-l, 
which follow from the periodic lateral boundary condi- 
tions, which require that e tkN = 1. The boundary con- 
ditions at infinity have a uniform gradient, i.e., 

lim ($(to + l,rc.) - $(to,?i)) = 1, n = 0, . . . , N - 1. 

m — >oc 

(1.9) 

Given the arbitrarily set of values <&(0,n), the solution 
for the row m = 1 is 



N-1 



*(l,n) = 1+ J2 $ (°> n ')3 



N(\n - n 



(1.10) 



n'=0 



where 



N-1 



9N( 



(") = n E e_K! C0S (M: n = 0, . . . ,N — 1, 

1=0 



(1.11) 

is the boundary Green's function, and hi corresponds to 
ki via the dispersion relation (1.8). The solution is given 
only for m = 1 , because we are only interested in the po- 
tential at sites near the interface. Note that the Green's 
function has the general property 



N-1 



E SN{n) = 1 



(1.12) 



n=0 



It is therefore good practice to check this normaliza- 
tion for each of the calculations presented below. Indeed, 
all our results obey this rule. 
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In general, the solution in the 'lower' regime is com- 
plicated by the variety of configurations. However, this 
solution is very simple for N = 2, when $(m, 0) is a lin- 
ear combination of e K,m and e _K/m . Since $(— j, 0) = 0, 
one is left with one unknown $(0,0), to be determined 
by the matching at row 0. 

For the special case N. = 2, the above procedure has 
led to the exact solution!] 



E. 



y{oo)e~ 



1,3 



-Eoo+l.oo I 1 



1- 



"f 3 



l+/3e 

-2k 



Of 3 







l+/3e" 



i>i, 



where 



oo+i.oo — hm Ej + i j 

j— too 



l+.g 2 (l)y(oo) 



< i < j - 2 
i=j-l 

i=j + l 
otherwise 

(1.13) 

0.5658..., 
(1.14) 



y(oo) = V3 - V2 = 0.3178..., e~ K f = 2- 
x/3 = 0.2679..., a = (1 + /3)g 2 (l)y(oo)/(2£ 00+ i, 00 ) = 
0.1281 ... and /3 = 5 - \/24 = 0.1010 . . .. For j = 0, the 
interface will transform into a step of size j = 1 with 
probability 1, hence E\$ — 1 and Eifl = for i ^ 1. The 
values of £!jj for < i, j < 4, up to the fourth decimal 
digit, are 



E 



0.4393 


0.5607 





0.3160 
0.1185 


0.5655 




0.3177 
0.0847 
0.0318 


0.5658 



0.3178 
0.0851 
0.0227 
0.0085 




(1.15) 



The first diagonal below the main, which represents the 
probabilities for the step to grow larger by one, Ej+ij, 
approaches its asymptotic value of -Eoo4 -i.no = 0.5658... 
exponentially, as the third row of Eq. ( 1.1 3| ) indicates. 
The diagonal above the main represents the probabili- 
ties for growths at the bottom of the fjord, Ej^ij, and 



corresponds to the second row in Eq. fll.l3p . These prob- 
abilities decay exponentially as the step size j grows. Ac- 
cording to the first row in Eq. (1.13), the elements Eij 
converge exponentially for large j's to a simple exponen- 
tial function: 



E. 



, yo = lim Eij = y(oo)e~ 



(1.16) 



These probabilities relate to the transition from a very 
large step to a step of size i. Next, the steady state vector 
P* is computed and used to evaluate the average upward 
growth probability (p up )*, which in turn, determines the 
average density p and the fractal dimension D. These 
computations are explained later in Sec. ||. 



Our previous paper does not specify details concerning 
the manner in which the system converges to the steady 
state in time. Besides addressing this issue, our present 
paper also treats DLA grown in wider geometrical peri- 
ods (still in cylindrical geometry). The basic approach is 
the same, i.e., we try to characterize the possible configu- 
rations of the interface for wider periods, and then write 
the evolution matrix, which is composed of the growth 
probabilities, which are computed from the Laplace po- 
tential, after proper normalization. The first difficulty 
encountered is in the characterization. For example, al- 
ready for a width of N = 3 one cannot characterize the 
interface using a single parameter as in the case N = 2, 
nor is it easy doing so using 2 parameters, or more. In- 
stead, we make a manual list of possible configurations of 
the interface, which we then order according to the dif- 
ference in height between the highest and lowest points 
on the interface. This difference is denoted by Am. Our 
order-O approximation includes only the configurations 
with Am < O. In our approximation, some of these 
configurations represent many other (excluded) configu- 
rations, in the sense that they have very similar growth 
probabilities, especially upward. This is because of the 
screening quality of the Laplace equation, which causes 
the potential to decay exponentially inside fjords. Thus, 
the deeper parts of the interface have a very small effect 
on the upward growth probability. The finite list of con- 
figurations is indexed arbitrarily, with an index usually 
denoted by i or j . Our experience shows that accurate re- 
sults are obtained, only when the order of approximation 
O is comparable to the width of the cylinder N. Thus, 
for wide periods, a high order calculation is called for. 
This causes the method to be ineffective for very wide 
periods, because the number of configurations grows ex- 
ponentially with the order of approximation. We con- 
ducted calculation up to N = 7. 

After selecting the finite list of configurations and 
obtaining the finite evolution matrix, we compute the 
steady state vector, which is the fixed point of the matrix 
(the normalized eigenvector with an eigenvalue of 1). For 
each configuration, we identify the upward growth pro- 
cesses (when the newly attached particle is higher than 
the rest). We then calculate the average upward growth 
probability (p U p)* as a weighted average over the config- 
urations. ^From (p U p)* we calculate the average density 
p and the fractal dimension D. The computed values of 
(Pup)* 5 from different orders of the approximation, are 
compared with numerical simulations in Table [|. 

In Sec. 

we introduce a simple Markov process, called 
the "frustrated climber" , which we solve exactly. A slight 
modification of the model is equivalent to site-DLA with 
a period of N = 2, which is presented in part II of this 
paperEI We then show a way of successively generaliz- 
ing the model to approximate bond-DLA with a period 
of N = 2 and with increasing orders O. We are able to 
check the approximations by comparing with the exact 
results of Rem. This model also enables us to investi- 
gate the rate of convergence to the steady state. In this 
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context we describe the convergence in terms of other 
eigenvectors, with eigenvalues whose absolute values are 
smaller than 1, and in terms of the infinite shift down op- 
erator. We show that the average upward growth proba- 
bility converges exponentially in time to its steady state 
value, with a characteristic time constant on the order of 
unity. In Sec. ID we generalize our method to cylindrical 
DLA with N > 2. We present in detail the calculations 
for N = 3 with = 1 and = 2, and for TV = 4 with 
= 1. Next we report on numerical results for wider 
periods and higher orders. In the final section we review 
the results and summarize. 



II. THE FRUSTRATED CLIMBER MODEL 

Consider someone trying to climb up a slippery infi- 
nite ladder. At each time step the climber climbs up 
one step with probability < p < 1, or falls all the way 
down with probability q = 1 — p. We call the climber 
"frustrated" , because the probability to get very high is 
exponentially small. We wish to compute the probability 
Pi(t) for the climber to be at height i after t time steps, 
for i = 0, . . . , oo. The Master equation for this problem 
is P(i + 1) = EP(t), where the matrix element Ei.j is 
the conditional probability that the climber moves from 
height j to i in a single time step. The rules of the model 
imply that 



i=j + l 

i = 
otherwise 



i>o, 



(2.1) 



so the matrix looks like this: 



E = 



q q q q 

p 
OpOO 
OOpO 



(2.2) 



This presentation helps us see the res emblance t o the 
dynamics of DLA with N = 2 in Eqs. ([Of, p..l5| ): Eq. 



(2.2) would approximate these equations if we were to re- 
place Ej + ij bypw £oc+i.oo and Eqj by q for all j, and 
neglect all other growth probabilities, which are indeed 
smaller. We shall discuss this and better approximations 
for DLA in the next subsections. Because the Markovian 
matrices for the two cases are similar for large j's, we ex- 
pect that some of the dynamical features are similar as 
well. We therefore present here an exact solution for the 
frustrated climber model, and then try to draw conclu- 
sions for generalized models which represent successive 
approximations for DLA. The advantage is that in the 
simple model of the frustrated climber it is possible to 
derive a simple analytic expression for the steady state 
and a complete description of the temporal convergence. 



The steady state equations for the frustrated climber 
model are 



P*i+i 
P 



P P*, i>0, 



(2.3) 



3=0 



-qp>, j>0. (2.4) 
One can easily check that this steady state is normalized, 



E^=E^= r 



= i. 



3=0 



3=0 



P 



(2.5) 



The average upward growth probability in the steady 
state is 

oo oo 

<fv>* = E p > p (j) = E p ip = p> ( 2 - 6 ) 

3=0 3=0 

where p up (j) stands for the probability to move upwards 
when the height of the climber is j. In this simple model 
Pup(j) =P for all j's. 

We now investigate the temporal convergence to the 
steady state. We define the vector v(i) by 



P(i) 



v(t). 



(2.7) 



Because P* and P(i) are probability vectors, Y^jLo Pj ~ 
S°°=o Pj(t) = 1, for any t, hence 



E^) = °- 

3=0 



(2.8) 



We substitute v into the dynamical equation and obtain 



P(* + l) =EP(t) = P* 
v(t + 1) = Ev(£). 



•Ev(f), 



(2.9) 
(2.10) 



Next, we look for the rest of the eigenvectors of the evolu- 
tion matrix (any eige nvector v with an eigenvalue A =/= 1 , 
has to obey Eq. (|2.8| )). Surprisingly, there are no eigen- 
vectors besides the steady state in this case. The eigen- 
vector equations are 



oo 

A«o = q E v 3 = °' 

3=0 

Awj+i = pvi(t), i > 0. 



(2.11) 



The first equation implies that either A = or vq = 0. 
In both cases, the last equation implies that v=0. 
We next introduce the infinite shift-down operator: 




10 
10 
10 



(2.12) 
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This operator causes a vector to "slide down" and inserts 
a zero at the evacuated component at the top. S has no 
eigenvectors at all, not even a fixed point (in spite of 
the fact that = 1 for j = 0, . . . , oo). In fact, 

Ev = pSv for all vectors v with J^JLo v j = 

Nevertheless, the convergence of P(t) to P* is simple. 
Starting from any initial state vector P(f = 0), the first 
application of E causes the first component to be set to 
its steady state value Po(t = 1) = q. At each subse- 
quent iteration another components is set permanently: 
Pi(t = 2) = qp, P 2 (t = 3) = qp 2 , etc. Pj becomes equal 
to P* after no more than j + 1 time steps. The context 
we are interested in is wider. We wish to compute the 
convergence of "observables" , i.e., the average of an arbi- 
trary function a(j) over configurations. We compute the 
average at time t 

oo oo 

(a)(t) ^aWPjit) = (ar + E#i(t), (2-13) 

3=0 3=0 

where (a)* = Y^Lo a (j)Pj IS tne steady state average. 
Starting from an initial deviation from the steady state 
v(0), each iteration causes a down shift and a multipli- 
cation by p, hence 

oo 

< a )(i) = (a)*+p t ^a(j+t)^(0). (2.14) 

3=0 



Equation (2.14) is the analogue of the standard eigen- 



vector description. We can also identify here the expo- 
nential decay of the factor p . For example, the function 
a(j) — Sjj g "measures" the probability of the climber to 
be at height jo (at any time). At time t the observed 
average probability is 



( a )(f)=i7 +pSo-t(o), 

for t<j , and (a)(t) = PT for t > j o 



A. First-order approximation for N = 2 



(2.15) 



We now return to Eq. (LIS), and try to approximate 
it by a sequence of models which are related to the frus- 
trated climber model. The simplest approximation would 
follow if we do not let the particle penetrate into th e fjor d 
at all. This is equivalent to setting Kf = oo in Eq. (1.13). 
According to these simplified rules, the particle can ei- 
ther stick at (0,0) and create a flat step of i = 0, or it 
can stick at (1,1) and increase the step height by 1. Let 
us denote the probability for the former event by q and 
the latter by p. In the first-order approximation we take 
p and q to be independent of the initial step size j, unless 
j = 0, in which case the step size increases with proba- 
bility 1. The Markovian matrix E for this case is almost 
identical to the case of the frustrated climber, 



E 



io q q q 

po o o o 

o p o o 

o o p o 



(2.16) 



the only difference being in the first column, where we 
denote qo — and po = 1. In part II of this paper we 
show that this model is exact for the case of site-sticking 
DLA for TV = 2LL3. 

The solution to this problem is very similar to that 
of the frustrated climber, with small modifications. The 
steady state is 



p; = p^pop>-\ j>i 



(2.17) 



where Pq can be determined using the normalization con- 
dition 

oo oo 

i = £p; = Po*(i+poX>')> 

3=0 3=0 
1-P 



=>Pn = 



(2.18) 



1 - P + Po 

The average upward growth probability is evaluated by 



(pVy = p *Po + (i-p *)p = 



Po 



1 - P + Po 



(2.19) 



The superscript (1) appears because it is the first-order 
approximation. We now need to choose p. One possible 
choice would be to take p = £00+1,00 = 0.5658, because 
this is the asymptotic upward growth probability, and 

then set q = 1 - p. This would give (jp^p)* rf 0.6973, 
to be compared with the exact value 0.68120. An al- 
ternative approximation would return to Eq. (1.13), but 
replace y(oo) by q, and then find q by solving 1 = p + q = 
[l+.92(l)g]/2+g. This yields p= 1-q = 2-^2 = 0.5858, 
and therefore (p$>* = V2/2 = 0.7071. 

We next calculate the average density and the fractal 
dimensionality—Similar to the argument used by Turke- 
vich and SchertHI, we consider a large number of growth 
processes n in the steady state. During this growth the 
aggregate would grow higher by h — (p up )*n. The to- 
tal volume covered by the new growth is hN d ~ 1 , where 
d = 2 is the Euclidean dimension. Thus, for N = 2 and 
for our first approximation the density is 



1 



P 



hN*- 1 (p up )*nN d - 1 (pup)*^- 1 



0.7171, 
(2.20) 



to be compared with the exact value p = 0.7340. Al- 
though our model does not really produce fractal struc- 
tures (due to the narrow width of our space), we can 
make an estimate of the fractal dimepsion in the same 
way Pietronero et al. estimated it inQiij. For a self sim- 
ilar fractal structure, one expects that a change of scale 
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by a factor N will change the average mass (number of 
occupied sites) of a N x N cut by a factor N D , where D 
is the fractal dimension. Assuming that the above proce- 
dure represents a coarse graining of the sites into N x N 
cells, we conclude that asymptotically 



p = N 



D-d 



(2.21) 



and this means that 

Hp) 



D = d 



ln(N) 



= l - ln j^p = 1.5202. (2.22) 



In Sec. IV we suggest a modified estimate of the frac- 
tal dimension, allowing for corrections to the asymptotic 
form ( pin] ). 

The study of the convergence to the steady state is 
again limited to the subspace of vectors v with X^jlo v o ~ 
0. The dynamic equation for z = is. 



B. Higher-order approximations for N — 2 



As mentioned earlier, the frustrated climber model re- 
in 



1.151) 



sembles the bond-DLA evolution matrix (1.13 
this section we approximate the full dynamics using in- 
creasingly more complex matrices. By doing so we do 
not improve on the accuracy of our previously published 
results^, but rather learn about the rate of convergence 
to the steady state. The method used in this section is 
generalized and applied to cylindrical DLA with wider 
periods in the next section. The case N = 2 is the sim- 
plest demonstration of this approach. 

The second-order approximation is to allow also tran- 
sitions of the kind j — ► 1 for j > 1. We also allow having 
arbitrary values in the top left 2x2 corner of the matrix , 
which we copy from the original matrix of Eq. (1.15), 



v (t+ 1) = q v (t) + ^2qvj(t) = (q - q)v (t), 

3 = 1 

v (t) = (q - «)*«b(0). (2.23) 

Since qa = 0, the exponentiated prefactor is negative, 
and therefore v (t) is oscillating during its decay. Af- 
ter the first iteration i>i(l) = pavo(0), regardless of its 
initial value. Afterwards it continues to follow vq, i.e., 
V\(t) — po(qo — q) t ~ 1 vo(0). After the second iteration 
^2(2) = popvo(0), and it also starts to decay exponen- 
tially with the factor (qo — q)- This happens for any 
j > 1; After more than j time steps (t > j) one has, 



v j (t)=p o p j - 1 (qo-q) t - j v o (0). 



(2.24) 



For short times and large indices t < j, the dynamics is 
governed by the shift down operator: 

00 

v(t)=v (Q)(qo-q) t h+p t Y / C 3 e U+t \ (2-25) 

3=1 

where e^' are the standard basis vectors, the components 
of the vector h are, 



h = 1, 

1 Po ( P 
hj = — 



p \qo-q 



(2.26) 



and the constants cj are determined by the initial condi- 
tions, 



c 3 = Vj(0) - v (0)hj, j = 1,2,. 



(2.27) 



For p > 0.5 the components of h explode exponentially. 
However, J2'jLo v j(^) ~ ® an< ^ therefore lim J -_ i . 00 Vj (0) = 
0. Thus, in order to cancel the divergence of the hj's, the 
Cj's must also explode exponentially and have an oppo- 
site sign. We note that because of this divergence h does 
not have a finite L\ norm and thus does not belong to 
the domain of E. Therefore it is not an eigenvector. 



E = 



9o 9i 9 9 9 

ro 7*1 r r r 

pi 

p 

p 



(2.28) 



We still require that the sum of the elements in each 
column be equal to 1, i.e., 



qo+r = 1, 
qi+n+pi = 1, 
q + r +p = 1. 



(2.29) 



In terms of standard DLA this means that we allow the 
particle to penetrate two sites into the fjord, but no more. 
Indeed it is exponentially improbable to penetrate deep 
into the fjord. This fact suggests a controlled approxi- 
mation for DLA. In each order of the approximation we 
allow the depth of penetration into the fjord to grow by 1. 
This is done by copyin g the (O + l)xO upper left block 
of the original matrix (1.13, 1.15), where O is the order of 
approximation. Asymptotic values are used outside this 
block, i.e., 

-^i+ij — -Eoo+1,00, j > O, 

Etj = y(oo)e- K f\ j > O, 1 < O - 2, 

n-2 



2/(00) 



i=0 
e -K f (n-l) 



3>0, 



(2.30) 



1 - e~ K f 

and the rest of the matrix elements are null. For exam- 



ple, in our case, O = 2, the constants in the matrix (2.26 
are 



G 



qo = o, 

ro = 1, 

6-3\/2 
q x = ^— = 0.4393, 

n = o, 

Pl = 3 ^~ 2 = 0.5607, 

q = y(oo) = V3 - \/2 = 0.3178, 
p = -Eoo+i.oo = 0.5658, 



r = y(oo) 



1 - e~ K / 



= 0.1163. 



First, the steady state is found by solving P* 
i.e., 

OO 

3=2 

OO 

p* = qoPo + qiPi + r Yl p h 

3=2 

pi=PiPi, 

/</'/• i>2. 

The solution to the last equation is 
P*=P&- 2 , j>2. 



(2.31) 
= EP*, 



(2.32) 



(2.33) 



Keeping this in mind it i s pos sible to exchange the two 
last equations of the set ( 2.32j ) with 



j2p*= p1 p?+ p j2p; 

3=2 3=2 



(2.34) 



Thus we obtain an autonomous and finite set of 3 
equations for 3 unknowns, namely, P *, P* and P 2 * = 
Sj°=2 Pj- The third parameter, P|, represents the total 
probability for the infinitely many configurations with 
j > 2. This reduction of the problem to three param- 
eters became possible because all of the configurations 
with j > 2 have exactly the same transition probabilities 
to the configurations j = and j = 1, and because they 
have exactly the same upward growth probability. Thus 
we obtain a fixed point equation for a 3 x 3 matrix, 



(2.35) 



It is guaranteed that a nontrivial solution exists, because 
the sum of the terms in each column of the finite ma- 



-p*- 




qo qi q 




-p*- 


PJ 




r n r 






. r 2 




pi p 




p* 



trix equals 1. Using the constants from Eqs. ( 2.31 ), the 
normalized solution obtained is. 



P 



<2) 
.(2) 



P^> = 0.2705, (0.2696), 



0.3184, (0.3113), 



P* (2) = 0.4111, (0.4191), 



(2.36) 



where the superscript denotes the order of approxima- 
tion and a comparison is drawn to the exact values in 
parentheses. By "exact" we refer to very high order cal- 
culations, or to values from simulations (which are the 
same up to the presented accuracy of 10~ 4 )l3. The ele- 
ments P* for j > 2 are evaluated using 



p; (2) = (i-p)p 2 * (2 V- 2 , j>2. 



(2.37) 



It is now possible to evaluate the average upward growth 
probability 



(P$>* = Po* r o + Pi> + P 2 *P = 0.6816, 



(2.38) 



where the exact valu e is .6812. The fractal dimension is 
evaluated as in Eq. ( 2.22j ), 



D 



(2) _ 



1.5530, 



(2.39) 



compared to the exact value 1.5538. 

The temporal convergence to the steady state in the 
second-order approximation can be analyzed using both 
the shift down operator and eigenvectors. The first eigen- 
vector of the matrix in Eq. (2.35) is the fixed point so- 
lution, which we denote by P*. Let us denote the other 
two (three-components) eigenvectors by h and g, and 
their corresponding eigenvalues by |Ao| > |AjJ. After t 
iterations of the evolution matrix we have 

P(t) = P*+c A h + ciA t 1 g, 



(2.40) 



where cq and c\ are constants determined by the initial 
conditions. The configurational average of some function 
a{j) with a(j) = <x(2) for j > 2, can be expressed in terms 
of these eigenvalues only, 



(a)(t) = (a)* + koXl + hXl, 



(2.41) 



where ko and fei are some other constants. A special 
function of this type is the upward growth probability, 
Pup(j) = (r ,pi,p,p,p,...). The eigenvalue with the 
largest absolute value other than 1, Aq, makes the largest 
contribution to the deviation from the steady state val- 
ues, and thus controls the temporal convergence. The 
characteristic time constant for the exponential conver- 
gence is, 



1 



(2.42) 



The eigenvalues obtained are Aq = - 
0.1257, using the constants of Eqs. 



-0.5599 and A^ ) 



( f2~3ll) . Hence, 



1.7. In order to describe the convergence of Pj(t) 



T <2) 

for j > 2 we use the vector v(t) = P(t) — P*, onc e mo re, 
and we perform a decomposition similar to Eq. ( [2.25| ): 



v(t) = coAoh + ciA^g 



3=2 



(2.43) 
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h 3 


= h jL 


9j 


= 9j, 


3 


= 0,1, 


hi 


= Pihi, 

/ \3- 2 


fj2 


= pm, 


3 


= 2, 


hj 




9j 


= . 


3 


> 2. 



where Co and Ci are the same as in Eq. (2.4C ) and the con- 
stants Cj for j > 2 are determined by the initial condition 
v(0). The vectors h and g are infinite generalizations of 
the finite vectors h and g, according to 



(2.44) 



Because p — -Eoo+i.oo > |Ao|, |Ai|, it is apparent that the 
components hj and gj diverge exponentially for large j's. 
This means that these vectors do not have a finite Li 
norm, and that they do not belong to the domain of E. 
Therefore, they are not eigenvectors, and A p and Ai are 
not eigenvalues of E. Nevertheless, Eq. ( |2.43| ) is still 
true. The effect of the shift down operator is manifested 
in the sum p* Cj^ +t - 

Using the same method it is possible to make higher 
order calculations. The steady state quantities resulting 
from the third order approximation are presented in Ta- 
ble EL in comparison with exact results. The eigenvalue 

— (3) 

with the largest absolute value is A = —0.5687, which 
has a greater absolute value than £?oo+i,oo = 0.5658. This 
means that a legitimate eigenvector exists for the infinite 
matrix. In the fourth and fifth order approximation we 

get A 4 ' 5 ^ ~ —0.5688. This suggests that the higher the 
order the more accurate is the evaluation of Ao and that 
the accuracy obtained is better than 10~ 4 . The typical 
time needed to settle in the steady state from any initial 
condition is therefore as short as 

t = 1.8. (2.45) 



III. DLA WITH N > 2 

The generalization of the exact methods from Ref.i to 
N > 2 is not straightforward. Trying to proceed along a 
similar line, one would try to parameterize the interface 
with a parameter i = 1,2,. . . , oo, and write the Master 
equation Pj(i + 1) = Yl'jLiEijPjit)- Unlike the case 
-/V = 2, the parameterization for N > 2 is very compli- 
cated. For instance, for the case N = 3 it is reasonable 
to try using two parameters, which indicate the height 
of two columns relative to the highest (or lowest) third 
column. However, this is insufficient because complex 
fjords (involving overhangs) might occur, as shown in 
Fig. |^. Instead of achieving a perfect para mete rization, 
we adopt the approximate approach of Sec. [IB, i.e., we 



take into account only a finite number of interface config- 
urations. These configurations are classified according to 
the maximum height difference between the highest and 
lowest particles on the interface Am. In the Oth-order 
approximation all the configurations with Am < O are 
included. The excluded configurations with Am > O are 



transformed into a configuration with Am — O, by fill- 
ing in the (O + l)th row below the highest particle; see 
Fig. [I| This transformation does not change the growth 
probabilities considerably. Especially, the upward growth 
probability would hardly change for large O. The vari- 
able Pi(t), where i corresponds to a configuration with 
Am = O, actually represents the sum of probabilities of 
all the configurations with Am > O, that have the same 
O uppermost rows, rather than represent the probability 
of the configuration i alone. This is analogous to P£ in 
the example above, see Sec. II B. After the finite set of 



configurations is chosen, the configurations are indexed 
with arbitrary consecutive numbers. Then, the growth 
probabilities for each configuration are computed by solv- 
ing the Laplace equation and by taking into account the 
bond multiplicity. Each growth process results in a dif- 
ferent final configuration, which must be identified with 
one of the configurations in the finite set. Special atten- 
tion is required for the upward growth processes, which 
might result in configurations with Am > O, which do 
not belong to the finite set. This is rectified by truncating 
the bottom row of the interface (considering it as fully 
occupied). The total upward probability for each con- 
figuration is added up and stored in a function p U p(i), 
later to be averaged over the steady state distribution of 
configurations. The growth probabilities are arranged in 
the evolution matrix, E, whose fixed point corresponds 
to the steady state distribution of configurations, which 
is required for evaluating (p U p)*, p and D. Because the 
matrix is finite, the existence of at least one fixed point 
is guaranteed. The other eigenvectors describe the rate 
of convergence to the steady state. 

The best way to demonstrate this approach is by show- 
ing a few sample calculations. The easiest ones are the 
first and second order approximation for N — 3 and the 
first order approximation for N = 4. After that we ex- 
plain the general algorithm for higher orders and widths, 
and report the results obtained numerically 



A. First order approximation for N = 3 

In the first order approximation we only look at the top 
row of the aggregate. For N = 3 there arc only 3 possible 
configurations (up to symmetry) , with the top row occu- 
pied by 1, 2 or 3 particles. Each configuration is indexed 
and for each configuration we identify the growth pro- 
cesses and the final configurations resulting from them; 
see Fig. |[ In part II of this paper we show that the 
calculation presented in this section can be used to solve 
exactly (no approximations) the case of site-DLA with 
N = 3. 

The first configuration (J — 1) grows upward with 
probability 1, thus _p U p(I) = I- The resulting configu- 
ration is i = 2, thus i?2,i = 1 and E^i = for i ^ 2. 
This concludes the construction of the first column of 
the evolution matrix. 
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In order to obtain the other growth probabilities we 
have to solve the relevant Laplace problems, for whic h 
we need the Green's function according to Eq. ( [l . 1 1 ) . 
For N = 3 we have h = ^-l for I = 0, 1,2. We recall 

that e~ K ( k ^ = q — \J q 2 — 1, where q = 2 — cos(fc)i and 
find that 



1, 



5- \/2T 



(3.1) 



The potentials of configuration j = 3 are described in 
terms of x = $(0, 1), as in Fig. 0. The Laplace equation 
is 



Ax = g 3 (0)x + 1 
6 - V2T 



= 0.2835. 



(3.7) 



There are 3 bonds leading to growth in site (1,0), which 
results in the configuration i = 1, hence 



and thus 



E 13 = ~x = 0.2835. 



(3.8) 



m o\ 1 - -93(0) V21-3 
33(1) = 53(2) = = . (3.2) 



These values obey the normalization condition ( 1.12 ). 

Because of the symmetry of the configuration j = 2, 
the potential can be expressed in terms of one variable 
x = $(0,0) = $(0,2), as shown in Fig. |. This kind 
of figure demonstrates the distribution of the potential 
$(m,n) over the lattice, and thus we call it a "po- 
tential diagram". The potentials $(1,0) = $(1,2) = 
1 + (1 — g 3 (l))x do not correspond to a growth pro- 
cess, but are important for solving for x. The potential 
$(1, 1) = 1 + 2xg 3 (l) corresponds to the upward growth 
process. The Laplace equation for x is 



4x = x + (1 - g 3 (l))x + 1, 
9 - V21 
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0.4417. 



(3.3) 



Growth in both sites (0, 0) and (0, 2) results in configu- 
ration i = 3, hence 



E 32 = A ~x= 18 - 2 f ri =0.5890, 
' 3 15 



(3.4) 



The upward growth process results in i = 2 after trunca- 
tion, and has probability 

Pu P (3) = E 2 . 3 = 1(1 + g 3 (l)x) = 0.7165. (3.9) 

The third element in the column is £3 3 = 0, which con- 
cludes the calculation of the elements of the evolution 
matrix, 



E (3,i) = 



0.2835 

1 0.4110 0.7165 
0.5890 



(3.10) 



where the superscript indicates that it is the first-order 
approximation for N = 3. The upward growth probabil- 
ities series is 



p up = (1,0.4110,0.7165), 



(3.11) 



which happens to be equal to the second row of the ma- 
trix. 

The normalized fixed point of the matrix is P* = 
0.0951, P 2 * = 0.5695 and P 3 * = 0.3354. The average 
upward growth probability is 

3 



5^i?JV(t) =0.5695. 



(3.12) 



where the numerator, 4, is inserted because there are 2 
bonds for each of the 2 growth sites, and the denomina- 
tor is the normalization factor N — 3. A growth process 
in site (1,1) results in an interface that does not belong 
to our finite set. In this approximation we only take into 
account the top most row of the interface, and therefore 
this interface is identified with configuration i — 2, i.e., 



We have performed some DLA simulations in the cylin- 
drical pgeometry for several values of N and measured 
(p U p) *tl The value obtained from simulations for = 3 
is 0.5462. The typical accuracy is on the order of 10~ 4 . 
The steady state average dens ity an d fra ctal dimension 
are evaluated using Eqs. ( |2~20| ) and fl2~22j ), 



3 15 y ' 

The transition to i = 1 is impossible, hence, E± t 2 = 0. 
It is easy to check that the second column of the ma- 
trix is normalized, i.e., Yli=i E%,2 = 1. The total upward 



growth probability for this configuration is 
Pu P (2) = £2,2 = 0.4110. 



(3.6) 



3<Pu P ) 



0.5853, (0.6103), 



D = 1 



ln((p up )*) 
ln(3) 



1.5125, (1.5506). 



(3.13) 



The values in parentheses are obtained from the same 
formulae, using the simulation value of (p up )*. The two 
other eigenvalues are c omplex, Ao,i = — 0.29 ± 0.28i, so 
according to Eq. fl2.42p t = 1.10. ' 
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B. Higher-order approximations for TV = 3 

The possible configurations of the interface in the 
second-order approximation are listed and indexed in Fig. 
^. The growth probabilities for the first 3 configurations 
were already computed in the previous section, but a re- 
arrangement of the upward growths is required in the 
evolution matrix. Now, the upward growth from config- 
uration j = 2 no longer stays at i — 2, but rather makes 
a transition to i = 4, and the upward growth from j = 3 
results in i = 5 instead of i = 2. Thus, we copy the pre- 
vious evolution matrix E^ 3 ' 1 ) into the upper left corner of 
the new matrix E^ 3,2 ) with the replacements: E%2 = 0, 



17,(3,2) 77.(3.1) 77,(3,2) „ j 77,(3,2 

hi a i — hj i , hi^ i — U, ana hi, 



E^ 1] . The un- 



4,2 — -^2,2 > -^2,3 ~~ u > ^5,3 

specified elements in the first three columns are all equal 
to zero. 

The next step is to go over each of the remaining con- 
figurations % = 4, . . . , 7, and compute their probabilities, 
which are inserted into the evolution matrix according 
to the final configuration in which the relevant growth 
process results. Configuration 4 is shown in Fig. 0. The 
Laplace equation is 

4y = y + x, 

4x = x + y + l + (l-g 3 (l))x, 



^x = —(7 - V21) = 0.5180, 
14 v 

y = x/3 = 0.1727. 

The growth probabilities are 
2 

E 6A = -x = 0.3453, 

4 4 
£5,4 = -!/=-^ 0.2302, 



(3.14) 



£4,4 — 



3" 9 
1 + 2xg 3 (l) 



= 0.4244. 



(3.15) 



The upward growth probability is p up (4) = £4,4 = 
0.4244. 

Configuration 5 is shown in Fig. [To| The Laplace 
equations are 

4y = y/4 + x + xg 3 (l) + yg 3 (0) + 1, 
4x = y + g 3 (0)x + g 3 (l)y + l, 

a- 

y = 0.4808, 
x = 0.4557. 



(3.16) 



The growth probabilities are 
E 7 . 5 = \x = 0.3038, 
£3,5 = I = 0.1603, 
£2,5 = |y/4 = 0.1202, 

E Ab = ±±9*m + V) = . 4157 . (3 . 17) 



The upward growth probability is p U p(5) = £4.5 = 
0.4157. 

Configuration 6 is shown in Fig. [ll| The Laplace 
equations are 

4y = y/4 + x, 
4.t = y + g 3 (0)x + 1, 

x = ^-(26 - 5a/2T) = 0.3067, 

!/=iit = 0.0818. (3.18) 

The growth probabilities are 





2 




-E-1,6 = 


r = 
2 


0.2044, 
8 


£3,6 = 


3 y = 


—x = 0.0545, 
45 


£7,6 = 


|y/4 


= x/15 = 0.0204, 






«? 3 (l)x) = 0.7206 


£5,6 = 





(3.19) 



E 



0.0 



The upward growth probability is p U p(6) 
0.7206. 

Configuration 7 is shown in Fig. [l2| The Laplace 
equations are 



Ax = x/4 + g 3 {0)x + 1, 
x = ^(21 - 4 V / 21) = 0.3051. 

The growth probabilities are 
2 

Ei, 7 = -1 = 0.2304, 

O 

£3,7 = \x/A = 0.0763, 

£5,7 = |(1+S3(l)a;)= 0.7203. 

The upward growth probability is p U p(7) 
0.7203. 

In summary, 



(3.20) 



(3.21) 



E5.7 



E (3,2) = 

0.2835 0.2044 0.2034 

1 0.1202 
0.5890 0.1603 0.0545 0.0763 
0.4110 0.4244 0.4157 
0.7165 0.2302 0.7206 0.7203 
0.3453 
0.3038 0.0204 



(3.22) 
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Pup = 

( 1, 0.4110, 0.7165, 0.4244, 0.4157, 0.7206, 0.7203 ) . 

(3.23) 

One can check that elements in each column of the ma- 
trix sum up to 1. Note that the majority of the elements 
are null. The normalized fixed point is, 

P* = ( 0.0685, 0.1011, 0.1145, 0.2680, 

0.2711, 0.0925, 0.0843 ) (3.24) 

with which we compute some steady state quantities, 



(Pup 



i=i 



0.5459, (0.5462), 



3(Pu P ) 



0.6106, (0.6103), 



D = 1 



ln ((Pup) 

ln(3) 



1.5510, (1.5506), 



(3.25) 



where once again, the values from simulation are shown 
in parentheses. It is apparent that the addition of con- 
figurations increases the accuracy of the results. The 
eigenvalues with the largest absolute values (except for 
1) are A ,i = -0.34 ± OAOi, hence r = 1.6. 

The third-order approximation yields 17 configura- 
tions. The final results are 

17 

<Pu P )* =^P;p up (i) = 0.5460, (0.5462), 



3(Pu F 



= 0.6104, (0.6103), 



D = 1 — 



ln((Pu P )*) 
ln(3) 



1.5507, (1.5506). 



(3.26) 



The eigenvalues with the largest absolute values (except 
for 1) are A ,i = -0.34 ± 0.40i, hence r = 1.6. 

It is interesting to inspect the histogram of the distri- 
bution of Pup{j), illustrated in Fig. [l^. One immediately 
observes that the upward growth probabilities are clus- 
tered in three groups: the top one at 1, the second just 
above 0.7 and the third, just above 0.4. It is easy to 
check that the top one corresponds to the configuration 
i = 1, the middle group corresponds to configurations 
that have two particles at the top row, and the bottom 
group corresponds to configurations with one particle at 
the top row. This suggests, that perhaps 17 different con- 
figurations are excessive, and the real number of effective 
configurations is around 3. An interesting question is 
whether it is possible to further reduce the number of 
configurations in higher-order approximations by includ- 
ing only "effective" ones. 



C. First-order approximation for N = 4 

Our last example is the case N — 4, for which we 
present the first-order calculation. First, we calculate 



the Green's function g 4 (n) according to Eq. (1.11). For 
N = 4, there are four possible values for k and k, namely, 



1, e 



k\ jy I 0, 2 J 7I "l2 7 '"' 6 

and e~ K2 = 3 — \/8. Hence, 

1 + 2(2 - V3) + 3 - V8 



V3, 



54(0) 
.94(1) 

.94(2) 



V3 + V2 
2 - = 0.4269, 



54(3) 



1 -3 + V8 V2-1 



4 2 
l-2(2-V3) + 3- V8_V3-V2 



0.2071, 



0.1589. 
(3.27) 



Once again, Eq. ( |l,12[ ) is obeyed. 

Figure [u] displays the relevant configurations. Con- 
figuration j — I grows into configuration i = 2 with 
probability 1, thus E2 1 = 1 and En — for i ^ 2. Also, 

Pup(l) = l- 

Configuration j = 2 is shown in Fig. [l5| The Laplace 
equations are 



Ax = y + g 4 (l)y + ( 9i (0) + g 4 (2)) x + 1, 
Ay = 2x + g 4 (0)y + 2g 4 (l)x + l, 

x = 0.5148, 
y = 0.6277. 



(3.28) 



The nonzero growth probabilities in the second column 
are £3,2 = fx = 0.5148, £4,2 = \y = 0.1569, and 
£2,2 = |(1 + 2 54 (1) x + 54 (2) y) = 0.3283 = Pup (2). 

Configuration j = 3 is presented in Fig. The 
Laplace equation is 



4x = x + (g 4 (0)+g 4 (l))x- 
>x = 0.4226. 



1- 



(3.29) 



The nonzero growth probabilities in the third col- 
umn are = = 0.4226 and £^,3 = 
I [1 + (34 (1) + 54 (2)) x] = 0.5774 = Pup (3). 

Configuration j = 4 is shown in Fig. [l7| The Laplace 
equation is 



ix = (g 4 (0) + g A (2)) x + 1, 
>x = 0.2929. (3.30) 



The nonzero growth probabilities in the fourth column 
are E bA = \x = 0.4393 and Ei A = | (1 + 2 54 (1) x) = 
0.5607 = p up (4). Note that this configuration already 
appeared for N — 2. 

The last configuration is shown in Fig. |l8|. The 
Laplace equation is 



11 



Ax = 1 = .94(0)2;, 
>x = 0.2799. 



(3.31) 



The nonzero growth probabilities in the fifth col- 



umn are E 



1,5 



0.2099 and E- 



2,5 — 



\ [3 + (34 (2) + 2 54 (1)) x] = 0.7901 = p up (5). This con- 
cludes the calculation of the 5x5 evolution matrix E^ 4,1 ) . 
The steady-state vector is 

P* = ( 0.0298, 0.4954, 0.2551, 0.0777, 0.1420 ). 

(3.32) 

It enables to calculate the following steady-state quanti- 
ties: 

<Pu P )* = P* 2 = 0.4954, (0.4657), 



4(Pu P ) 



0.5046, (0.5368) 



D = 1 



ln((pu F 



ln(4) 



1.5066, (1.5512), 



(3.33) 



where again, the values in parentheses are from simu- 
lation. The eigenvalues with the largest absolute value 
after 1 are A0.1 = — 0.16± 0.38i, hence r = 1.1. 

It is also possible to conduct these calculations using 
different boundary conditions at the bottom; rather than 
assuming that there is a filled row of occupied sites below 
the configuration, it is possible to assume that each unoc- 
cupied site at the lowest row of the configuration is above 
an infinite fjord that extends all the way below. The two 
possibilities are explained in Fig. [l^. Performing the 
calculations with infinite fjords is a bit simpler, because 
there are less configurations, e.g., the configuration i = 4 
would-jjot appear in the first-order approximation for 
N = M. 



D. Higher order computations 



As one increases N and the order of approximation 
O, the number of configurations increases exponentially, 
and it becomes harder to go over all of them manually. 
However, it is possible to construct a computer algorithm 
to perform the procedure described here. The main chal- 
lenges are the automatic configuration recognition and 
automatic computation of the exact growth probabilities 
per configuration. In this section we explain the algo- 
rithm and report some of the important results. 

The algorithm follows the method outlined in the ex- 
amples of the previous sections, i.e., it goes over all the 
possible configurations of the interface. In the sample 
calculations we have initially made a list of all the possi- 
ble configurations, called the index. Instead of doing this, 
the program starts with only one configuration, namely 
the flat one (all the sites of the top row of the aggregate 
are occupied), which is indexed by j = 1. This configura- 
tion grows with probability 1 to a new configuration that 



has one particle at the top row, while the row below it 
is fully occupied. This new configuration is inserted into 
the list of configurations with an index j = 2. Therefore, 
the program sets i?2,i = 1 and p U p(l) = 1- Then the 
program continues by handling the next configuration in 
the list, namely j = 2. For each configuration, it solves 
the Laplace equations and calculates the growth prob- 
abilities. Each growth process may create a new con- 
figuration. The resulting configuration is first checked 
for consistency with the desired order O; configurations 
which have Am > O are truncated, as in Fig. ^. One 
then compares each 'new' configuration with the existing 
list of configurations. If it does not exist in that list it is 
added at the end of the list, and indexed consecutively. 
If the index of the configuration that results from the 
growth process is i and the index of the initial configu- 
ration is j then the growth probability is inserted into 
the matrix element Eij. The total sum of all the up- 
ward growth probabilities of the initial configuration j is 
stored in p U p{j)- The main loop stops when the program 
finishes to process the last configuration in the index list. 
At this stage the Markovian evolution matrix E is irre- 
ducible and closed, i.e., J2i — 1 f° r every j. Then the 
fixed point P* is calculated, by taking an initial vector 
and iterating E on it many times until it converges (for 
very large matrices this is much faster than using any of 
the MATLAB library functions). The average upward 
growth probability is calculated using 



(p up )* = ]Ti> up (j), 



(3.34) 



the average density and the fractal dimension are then 
comp uted using the left hand side of Eq. (2.20) and Eq. 
( PI )- 

One of the challenges of the computer algorithm is the 
recognition of configurations. This recognition is impor- 
tant so that each growth process will be inserted into the 
evolution matrix Eij with the correct index i (j is the in- 
dex of the configuration before growth). The recognition 
maybe difficult because configurations that seem differ- 
ent may actually be equivalent. By equivalent we mean 
that they have the exact same set of transition (growth) 
probabilities. The solution to the Laplace equations is 
determined uniquely by the shape of the interface, there- 
fore all of the configurations with the same external in- 
terface are equivalent. The description of the interface is 
not a trivial task though. We find that an efficient way 
to characterize an interface is by the set of empty sites 
that are connected to infinity. Of course, it is sufficient 
to specify only empty sites that are not higher than the 
highest particle in the aggregate, because all of the empty 
sites above it are connected to infinity. Figure [2(] shows 
an example of two configurations that are not identical, 
but they have the same exterior contour. Both of them 
have a single empty site that is connected to infinity. 

In order to reduce greatly the number of configura- 
tions it is advisable to take symmetry into account, i.e., 
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all the configurations which can be obtained from one 
another using a rotation around the axis of the cylinder 
have the same growth probabilities and the same steady 
state weights. The same is true for mirror images. In- 
stead of taking all of them into account, we choose one as 
a canonical representative of the whole set of symmetric 
configurations. 

The results are summarized in Table [|. By comparing 
the approximations to accurate results from simulations, 
it seems that in order to obtain a relative accuracy of 
about 10 -3 one has to use at least an order of approx- 
imation of O = N — 2 (except for N — 3, where one 
still has to use the second-order approximation). This 
becomes very difficult already for N = 6, where in the 
fourth-order calculation there are 49678 different config- 
urations up to symmetry. 



IV. DISCUSSION 

This paper treats DLA as a Markov process. The 
Markov states are the possible shapes of the interface, 
and the Markovian evolution matrix E is calculated an- 
alytically using exact solutions of the Laplace equations, 
with proper normalizations. We propose a truncation 
scheme that takes into account only a finite number of 
states. The states are ordered according to the maxi- 
mal difference in height between the highest and lowest 
points on the interface, Am, and in each order of trun- 
cation O, only the states with Am < O are included. 
We justify this approach by the fact that the potential $ 
decays exponentially in deep fjords, and thus the shape 
of the interface in its deeper parts has very little effect 
on the growth probabilities. We perform this calculation 
for N — 2, and verify that indeed it converges to the 
known analytic solution. We adopt the same approach 
for higher values of the width N, between 3 and 7, and 
calculate the average density p in good agreement with 
simulations. The fact that the number of configurations 
grows exponentially with N and with O, makes the com- 
putation less effective than simulation for large N. 

We observe that the method converges as a function 
of O, also for higher values of N. Let us denote the cal- 
culated average steady-state density of an aggregate of 
width N in the O'th-order approximation by p c (N,0). 
We observe that p c (N, O) converges to a finite limit very 
rapidly as a function of O. In fact, a relative accuracy of 
10" 3 is achieved for O = N - 2 (except for N = 3). This 
enables us to obtain accurate results for 3 < N < 6. The 
drawback of this method is that the number of configura- 
tions diverges exponentially with O and N, and therefore 
it is possible to perform the calculations only for rela- 
tively low TV's and O's. Our computer was strong enough 
to perform the calculation only in the third-order approx- 
imation for N — 7, and therefore the result for N = 7 is 
not very accurate. One would hope that it may be pos- 
sible to perform low-order approximations for large iV's 



and then extrapolate, in order to estimate the results for 
large O's. Indeed, it is reasonable to conjecture the scal- 
ing law p c (N,0) = p(N)f{N/0), where p(N) is the ex- 
act (O — > oo ) density, as a function of N, and f(N/0) is 
a universal scaling function that obeys ]im x —tof(x) = 1. 
Our investigation shows that in spite of the fact that the 
conjecture is not very accurate for O = 1 and O = 2, it 
is quite good for higher values of O, and presumably also 
for higher values of N. This scaling relation may help to 
perform the extrapolation O — > oo for higher values of 
N. Paradoxically, it is very hard to obtain data points 
for large TV's and O's, and thus to extract the scaling 
function accurately. Thus we are unable to make the ex- 
trapolation even for N — 7, and we estimate p(N) by 
the highest-order approximation available. However, we 
suggest an alternative way to obtain p c (0,N), namely 
by simulation: it is possible to perform a regular DLA 
simulation in cylindrical geometry, only that one has to 
keep the O'th row below the highest particle in the ag- 
gregate constantly filled. Measuring the average density 
of the aggregate in such a simulation would approximate 
p c (N,0). This simulation would be faster than a regu- 
lar simulation, because particles would stick faster, due 
to the fact that they have less free space. This study 
would perhaps yield the scaling function f(N/0), and 
enable extrapolation of lower order approximations for 
higher iV's, should anyone venture to perform them on 
more powerful computers. In light of this discussion we 
suggest a more efficient way to perform DLA simulations 
in cylindrical geometry. We argue that one can obtain a 
relative accuracy of 10 -3 if one follows just the N—2 top 
most rows of the aggregate. This should save some time, 
because the diffusing particle would stick faster, and it 
would also require less memory. This is not to say that it 
is sufficient to grow the aggregate until it reaches a height 
of N — 2, but rather, to perform many more growth pro- 
cesses, and each time the aggregate reaches a height of 
N — 1, truncate the bottom row. 

We also discuss the temporal rate of convergence of the 
system to its steady state. In this context we find that 
there is an exponential convergence to the steady state, 
and we calculate the characteristic time constant r. This 
is demonstrated using the simple model of the frustrated 
climber. The convergence is described in terms of the 
eigenvalues of the Markovian matrix, and in terms of the 
infinite shift-down operator. 

Considering the fractal dimension, Pietronero et al. 
sugge sted that p(N) = N D ~ d , as mentioned in Eq. 



( [2.21 ). In principle, one should always include an am- 



plitude and finite size corrections of the form 
p(N) = AN- a (1 + B/N + ...) 



(4.1) 



where a = d — D, and A and B are constants. The sec- 



ond term appearing in Eq. (4.1) is a correction to scaling 
term. Generally, there is an infinite sum of such terms 
with higher negative powers of N. Because we have data 
only for small values of N, these correction terms may 
be large, but since we have only a few accurate data 
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points (p(N) for N — 2, 3, . . . , 6), we try to extract the 
parameters a, A and B only, and not higher order terms. 
Using the three results for N = 4,5, 6, we determine the 
three unkown parameters to be A = 0.82, B = 0.35 and 
a = 0.362, hence D — 1.64. The deviation from the well 
know value of D = 1.66 can be attributed to system- 
atic error due to the omission of higher order finite size 
correction terms. We fit simulation datao for N = 3, 
4, 5, 6, 7, 32, 48, 64, 96, 128, to a higher-order approxi- 
mation p(N) = AN~ a (1 + B/N + C/N 2 ) , and find that 
C = -0.205, B = 0.561, A = 0.761 and a = 0.339, which 
means that D = 1.661. The maximum relative error of 
the fit is 1.2 x 10~ 3 , and the average relative error is 
1.0 x 10~ 3 , which is in good agreement with estimated 
accuracy of the simulations. 
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FIG. 1. The coordinates (m, n) describe the location on a lattice that is two sites wide. The gray sites belong to the interface 
of the aggregate, which is shaped as a step of size j. 



Upward growth: 
1 bond 




FIG. 2. Possible growth processes that change the interface from an initial step size j — 3 to a final size i = 4,0, 1, 2. The 
growth probability is determined by the potential and the number of bonds associated with the site where growth is to occur. 
Ei j is the conditional probability to grow from an initial step size j to a final step size i. The normalization follows from Eq. 
(1.10). 
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FIG. 3. An example of an interface configuration for N = 3 that cannot be characterized using the height differences of two 
columns relative to the third. 
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FIG. 4. Configuration (a), with Am = 2, is truncated by taking only the top row, and turns into configuration (b), with 
Am = 1, in the first-order approximation (O = 1). 
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FIG. 5. The three possible configurations in the first-order approximation for N — 3, up to translation symmetry. The 
arrows indicate the possible transitions due to growth processes. The transition probability from configuration j to i is denoted 
by Ei.j . 



1+2g(1)x + 

3 



m=1 
m=0 



j=2 




1+(1-g(1))x 

3 



m 

n= 0, 1, 2 

FIG. 6. A "potential diagram": the potentials "I>(m,n) of the configuration j = 2, expressed in terms of the variable x. 
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FIG. 9. The "potential diagram" for configuration j = 4. 
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FIG. fO. The "potential diagram" for configuration j = 5. 
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FIG. 11. The "potential diagram" for configuration 
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FIG. 12. The "potential diagram" for configuration 
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FIG. 13. The distribution of p up over configurations for the third-order approximation for N = 3. 
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FIG. 14. Possible configurations in the first-order approximation for N 
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FIG. 16. The "potential diagram" for configuration 
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FIG. 17. The "potential diagram" for configuration 
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FIG. 18. The "potential diagram" for configuration j = 5. 
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FIG. 19. The two top rows of a configuration are shown in (a). Two possible extensions for the rest of the configuration 
below are (b) , with a filled row right below the configuration (this boundary condition is used in the calculations presented in 
this paper), or (c), with the bottom row of the configuration repeating itself ad infinitum, creating an infinite fjord. 
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(a) (b) 



i 




FIG. 20. Even though configuration (a) and (b) are not identical, they are equivalent because they have the same growth 
probabilities. Both configurations have the same external interface contour, which is characterized by the set of sites that are 
connected to infinity. In this example there is only one such site, which is not higher that the aggregate, and it is marked by a 
circle 



TABLE I. The two-dimensional approximate results for various channel widths N and for different orders of approximation 
O. The quantities presented in each table cell are the average upward growth probability (p up )* and the number of configurations 
N c . The approximate results are compared with simulations. 



N/O 


simulation 


1 


2 


3 


4 


5 


6 


3 


0.5462 


0.569489 
3 


0.545911 
7 


0.546046 
17 


0.546126 
45 


0.546132 
127 


0.546132 
371 


4 


0.4657 


0.495435 
5 


0.464571 
20 


0.465395 
98 


0.465730 
575 


0.465765 
3640 


0.465768 
23676 


5 


0.4106 


0.444088 
7 


0.407582 
47 


0.409497 
457 


0.410414 
5539 


0.410547 
69791 




6 


0.3696 


0.405619 
12 


0.364352 
131 


0.367295 
2217 


0.369172 
49678 






7 


0.3377 


0.375448 
17 


0.330112 
337 


0.333622 
10403 









TABLE II. Some steady state results of the third order approximation 









p. 


Pi 


P* 


PI 


r; 


3rd order 


0.6812 


0.2696 


0.3114 


0.1820 


0.1029 


0.0582 


0.0329 


Accurate 


0.6812 


0.2696 


0.3113 


0.1809 


0.1032 


0.0586 


0.0332 
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